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ABSTRACT 

We present recent results of 3D magnetohydrodynamic simulations of neutron 
stars with small misalignment angles, as regards the features in lightcurves produced 
by regular movements of the hot spots during accretion onto the star. In particular, we 
show that the variation of position of the hot spot created by the infalling matter, as 
observed in 3D simulations, can produce high frequency Quasi Periodic Oscillations 
with frequencies associated with the inner zone of the disk. Previously reported sim- 
ulations showed that the usual assumption of a fixed hot spot near the polar region 
is valid only for misalignment angles 6 relatively large. Otherwise, two phenomena 
challenge the assumption: one is the presence of Rayleigh-Taylor instabilities at the 
disk-magnetospheric boundary, which produce tongues of accreting matter that can 
reach the star almost anywhere between the equator and the polar region; the other 
one is the motion of the hot spot around the magnetic pole during stable accretion. In 
this paper we start by showing that both phenomena are capable of producing short- 
term oscillations in the lightcurves. We then use Monte Carlo techniques to produce 
model lightcurves based on the features of the movements observed, and we show 
that the main features of kHz QPOs can be reproduced. Finally, we show the behavior 
of the frequencies of the moving spots as the mass accretion rate changes, and propose 
a mechanism for the production of double QPO peaks. 

Keywords: accretion, accretion discs; instabilities; MHD; stars: neutron; stars: oscil- 
lations; stars: magnetic fields 



1 INTRODUCTION 

Quasi Periodic Oscillations (QPO) have been observed in 
many X-ray sources for the last 24 years ( van der Klis et al."| 
1985| |Lamb et al.|1985| ). They appear as broad peaks in the 
power spectra, with a behavior dependent on the source state, 
the mass accretion rate, and other physical properties of the 
sources. They are present in both neutron star (hereafter NS) 
and black hole (BH) sources, sometimes with peculiar simi- 
larities between the two classes of objects (see |van der Klis| 
|2006| for a review). 

Of particular interest are the kHz QPOs: they were first 
discovered in the Low Mass X-ray Binary (LMXB) and Z source 
Sco-Xl ( [van der Klis et al.[T 996 ), to be then observed in the 
spectra of nearly all Z and atoll sources, as well as several 
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weak LMXBs (see for example Strohmay er et al.||l996[ |Wi-| 
jnan ds et al.||1997l). They hav e frequencies between 300Hz 
and 1.2kHz (see |van der K lis 2006 ). The peaks move with 
timescales of hours/days and their quality factor Q = v/Av 
(ratio between the frequency and the width of the peak, a 
measure of coherence) is not constant, but can be related with 
the QPO frequency ( [Mendez et al.|200H|di Salvo et al.|2003t 
Barr et et al.|2006 ) and with the mass accretion rate (see for 
example Mend ez]2006| ) . 

The kHz QPOs usually appear as a pair of peaks (the 
twin peaks), labeled "upper" and "lower" (and their frequen- 
cies, accordingly, v u and v,). While on timescales of hours the 
two frequencies change, their difference Av = v u - v l tends 
to be almost costant. As with the discovery of burst oscilla- 
tions ( Strohm ayer et al.||1996|) and Accreting Millisecond X- 
ray Pulsars (hereafter AMXP) ( Wijnands & van der Klis|1 998 ) 
the spin frequencies of more and more accreting neutron 
stars were measured, an interesting behavior was observed: 
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Av was usually found to be close to either v* or 0.5 v*. Recent 
works (see for example Mendez & Belloni 2007 ) call this as- 
sumption into question. In particular, it can be shown that Av 
could in principle be independent from v*. 

It has been shown that the kHz QPOs frequency is cor- 
related, at least on a given system and for short timescales, 
with the X-ray count rate (Mend eFet al.|19 99 ) and thus also, 
probably, with the mass accretion rate. 

Many models have been developed to explain the kHz 
QPO phenomenon. They involve a number of different mech- 
anisms, like for example relativistic precession frequencies 
([Stella & Vietri 1998), relativistic resonances ( |Kluzniak~&[ 
|Abramowicz||200i' ), beats between the keplerian frequency 
at a particular radius in the disk and the frequency of the 
star ||Alpar & Shaham|| 1 985 \ [Lamb et al.|[l985] |Miller et al.| 



One of the interesting phenomena observed in 3D 
simulations is the movement of spots around the magnetic 
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Lamb & Miller 2004 ), oscillations in the comptonizing 



medium (Lee & Miller 1998). Recently, Lovelace & Romanova 
( |2007| ) and Lovelace et al. ( |2009| ) found and analyzed a ra- 
dially localized Rossby wave instability which exists in the 
inner region of a disk around a magnetized star where the 
angular rotation rate of the matter decreases approaching the 
star. This instability may explain the twin kilo-Hertz QPOs ob- 
served in some LMXBs. The importance for QPOs of the disk 
angular rotation rate having a maximum was discussed ear- 
lier by |Alpar & Psaltis[ ( [2005] [2008] ) 

Most of the models of kHz QPOs, based on the fact that 
QPOs also appear in black hole systems and there seems to 
be a link between v { with frequencies observed in black hole 
systems (Psal tis et al.| [l999 ), try to explain the behavior of 
kHz QPOs in both the neutron star and the black hole cases, 
excluding the possibility of a surface emission. Gilfanov & 
Revnivtsev] ( |2005| ) claim that, according to spectral emission 
properties of QPOs, in Atoll and Z sources periodic and quasi- 
periodic emission is more likely to be produced on the surface 
of the star. 

During the past years MHD simulations have been 
amongst the most fruitful methods to investigate accreting 
magnetized stars (see for exam ple [Goodson et al.|1997HRo-| 



manova et al. 2002 , 2003 , 2004; Kulkarni & Romanova 2005 ; 



Bessolaz et al.|2008|[Long et al.|2 008 ; Ro manova et al.|2 008 ). 



3D simulations have given insight into the way matter ac- 
cretes, showing in more detail trajectories, shapes of magnetic 
field lines, shapes and movements of hot spots on the surface 
of the star (see for example Romanov a et al.|2003[|2004] ). 

Simulations show that matter may accrete in several dif- 
ferent regimes: 

• the stable accretion regime, when matter is stopped by 
the magnetic field and conveyed to the magnetic poles via the 
so called funnel flo w ( [Romanova et al.|2002[ |2003| |Kulkarni| 
|&Romanova|2005| ); 

• the unstable regime, when Rayleigh-Taylor instabilities 
form at the boundary between the disk and the magne- 
tosphere (Roman ova & Lovelace | [2006} |Romanova et al.| 
2008 ; Kulkarni & Romanova 2008}, and accretion takes place 
through tongues of matter from the inner edge of the disk to 
random places between the equator and the magnetic poles; 

• the magnetic boundary layer regime, during which the 
disk is very close to the surface of the star and the accretion 
takes place via two big antipodal equatorial tongues which ro- 
tate with the frequency of the inner disk ( Rom anova & K ulka- 
|rni|20091 ). 



pole of the accreting star, both in stable (Romanova et al. 
200 3] 120041 |2006) and in unstable regime ( [Kulkarni & 
Romanova 2009|). It was noticed that in some cases, such 
as in the magnetic boundary layer regime, the moving spots 
produce QPO features ( Romanova & Kulkarni 2009 ). Usually, 
models of NS surface emission (like in the case of Accreting 
Millisecond Pulsars) consider funneling of matter onto the 
polar caps as a static process, with a fixed polar cap and thus 
a fixed emission spot 

citep[e.g.][]Beloborodov:2002p4173,Poutanen:2003p4169. 
3D simulations show that this may not be the case in 
particular for small misalignment angles. 

Recent works ( La mb et al.|2009a]b] ) investigate the wan- 
dering of hot spots at which matter accretes in sources with a 
small misalignment angle, finding that it might affect the co- 
herence of the pulsed signal from neutron stars and even de- 
stroy it, being the origin of the transient emission in AMXPs. 
The movements hypothesized in those works, though, are 
different from the movements observed in our 3D simula- 
tions, where the motion of the hotspots is a rotational mo- 
tion around the magnetic pole or in the equatorial region, 
on timescales of milliseconds. Moreover, the dimension of 
the hotspot in those two works (a radius of ~ 25°) does not 
find confirmation in our simulations. It corresponds to the di- 
mension of the whole "hot ring" where our hotspots (much 
smaller) move (see Fig. [TJ . 

In this paper we show that the observed movements carry 
with them the information of disk frequencies near the inner 
disk, and their effect is the production of quasi-periodic os- 
cillations at the relative frequencies. We perform a series of 
3D MHD simulations (mostly for the cases of very small 0) 
and analyze the rotation of the spots and the related frequen- 
cies in detail. MHD simulations can only follow the motion of 
the disk on very short timescales compared to the disk evolu- 
tion timescale and the observational time. For this reason, we 
use the features of the moving spots produced in simulations 
to generate long lightcurves by means of Monte Carlo tech- 
niques, showing that the motion observed produces QPOs. 
Then, we show that the phenomenon follows a predictable 
behavior as the mass accretion rate changes. 

We start in ^2]with a general description of the 3D MHD 
model used for this work. Then, we show the way spots move 
in our simulations in ^3) We then show with a Monte Carlo- 
generated lightcurve our model of surface QPO production in 
^4) Finally, in ^5] we show how the angular velocities of hot 
spots correlate with the mass accretion rate. 



2 MODEL USED FOR 3D SIMULATIONS 

We model the accretion disk around a neutron star as a sin- 
gle fluid plasma, interacting with the star's magnetic field. We 
then solve the complete system of MHD equations in the co- 
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ordinate system rotating with the star (Kol doba et al.|2002] ) 

(1) 



_^+V(pv) = 



g(pv) 

at 



+ V-^ = pg + 2(pv)xfi-pfix(fixR) (2) 



Jt 



V x (v x B) 



dt 



■ V • (pSv) = o, 



(3) 
(4) 



where p is the density, v is the velocity of the fluid, g is the 
gravitational acceleration, B is the magnetic field, St is the 
momentum flux-density tensor, and S is the density of en- 
tropy. The equations are solved by means of the Godunov 
scheme (Kulikovskii et al. 200T]), using a cubed-sphere grid. 
The outward transfer of angular momentum in the disk is pro- 
duced by an a-type viscosity (Shaku ra & Sunyaev||1973| ). A 
detailed description of these equations and their implementa- 
tion can be found in |Koldoba et al.| ( [2002l ). 

The potential well around the neutron star is represented 
by the pseudo-Newtonian potential described in Paczyhski & 
Wiita]( [T980] ), in order to take into account relativistic correc- 
tions to the movement of matter around the star. This poten- 
tial is described by the equation 



V(r) = - 



GM 
r-R a 



(5) 



where R g = 2GM/c 2 is Gravitational Radius (equal to the 
Schwarzschild Radius). This form of the potential gives the 
correct values for the marginally bound circular orbit and in 
particular the innermost stable circular orbit (ISCO) in the 
Schwarzschild non-rotating BH case. This can be useful to in- 
terpret the results of this work in terms of theoretical models 
taking the ISCO as a fundamental quantity for understanding 
QPOs. 



2.1 Conventions and units 

All quantities used in the 3D code are dimensionless. Here- 
after, given a dimensional quantity ^, we call 9/ = tne 
corresponding dimensionless unit, where is the reference 
value for . Every result can be applied to different physical 
systems, evaluating the %s by choosing the right values of 
three fundamental dimensional parameters: 

• M: the mass of the star; 

• R: its radius; 

• its magnetic moment, corresponding to a dipolar sur- 
face magnetic field B = (jl/R 3 . 

At the start of simulations, we fix the following dimensionless 
parameters: 

• the corotation radius r c , defined as the radius where the 
Keplerian angular velocity equals the star's angular velocity. It 
sets the rotational velocity of the star (s ee Eq.[6)); 

the viscosity a, on the model of [Shakura & Sunyaev 



( [1973] ) 

• the magnetic moment Jl, and possible non-dipolar com- 
ponents, which are not used in this work. Alternatively, Jl can 
be considered as a parameter to set the size of the magneto- 
sphere and the mass accretion rate (see Eq.[8]); 



• the misalignment angle 0, the angle between the mag- 
netic axis and the rotation axis. 

The reference quantities are calculated as follows: 

• M = M; 

• r = .R/0.35, the reference distance from the center of 
the star, chosen this way for convenience, because of the vicin- 
ity of r = 1 to the usual range of magnetospheric radii; 

• jU = ju/jCt and B = jU /r 3 ; then, the surface magnetic 
field is = (0.35)- 3 B oJ u; 

• v = (GM/r ) 1/2 , the reference velocity; 

• co = (GM/r 3 ) 1/2 , the reference angular velocity; 

• P = 2nr /v , the Keplerian period at r = 1 as the refer- 
ence time; 

• p = B l/ v o as tne reference density; 

• M = p v r 2 as the reference mass accretion rate; 

• L = M Vq is the reference luminosity 

from which we obtain 

(0.35) 3/2 (GM) 1/2 



2tt (r c R) 3 / 2 



1 



(0.35) 7/2 fB, 



M = M n M = - f ^ ) R 5/2 M 

(GM) 1/2 V U 



(6) 



(7) 



(8) 



35) 9/2 (GM) 1/2 ^j 2 R 3/2 L (9) 

where v*, t*, M and L are the frequency of the star, its period, 
the mass accretion rate and luminosity respectively. M and L 
are returned by the code at every time step. L is calculated un- 
der the assumption that all the energy of the accreted matter 
is converted into luminosity. Equation [8] gives us a way to use 
the parameter jl to effectively measure the accretion rate. In- 
stead of using the simplest interpretation, according to which 
p is a dimensionless measure of the magnetic field, we can in 
fact change B in order to make B* constant, changing effec- 
tively the mass accretion rate reference value M . 



2.2 General information about the simulations in this 
paper 

Unless otherwise specified, throughout this paper all results 
are referred to a neutron star with R = 10km, M = 1.4M , 
6 = 2°. The magnetic field will be either B = 10 8 G or B = 
5- 10 8 G. With these parameters, the gravitational radius of the 
Paczynski-Wiita potential is R g w 4.16km while the ISCO is at 
a radius R IS co w 12.5km (in dimensionless units, r = 0.435). 

The grid used for most of the simulations has 72 cells 
in radial direction and 31 along every side of the six blocks 
of the cubed sphere. The simulations shown have been pro- 
duced in the period January - June 2009, using the High Per- 
formance Computing facilities of NASA (Pleiades, Columbia) 
and CINECA(BCX). 



3 3D SIMULATIONS AND SPOT MOVEMENTS 

In this section we show how spot motions can produce visible 
features in the light curves. 
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Figure 1. Evolution of the position of the hot spot during accretion, for the case r c = 1.5 (t* = 4.1ms), a = 0.04, \x = 2, a good example of stable 
accretion. The time difference between the files is 0.5P . In this case all accretion takes place around the magnetic pole. 
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Figure 2. Evolution of the position of the hot spot during accretion, for the case r c = 1.8 (t* = 5.4ms), a = 0.03, fx = 0.9, example of unstable 
accretion. The time difference between the files is 0.25P . In this case matter pressure is such that the magnetic field is not able to convey all 
accreting matter to the magnetic poles. Instead, instabilities form at the truncation radius, creating tongues of accreting matter which fall onto 
the surface in random places. Accretion takes place both in the polar region and in the rest of the surface 



3.1 Stable, unstable, strongly unstable regimes 

The magnetic field of a NS is usually strong enough to trun- 
cate the accretion disk at a certain distance, called inner ra- 
dius or magnetospheric radius, where the magnetic and mate- 
rial stresses are equal ( Ghosh & Lamb |1979| ). In most cases, if 
the inner radius is comparable with the corotation radius, ac- 
cretion takes place in a stable fashion: matter is stopped in the 
disk by the magnetic field and moves towards the magnetic 
poles with the magnetic field lines acting like an invisible fun- 
nel (see Fig. [I] left) . Depending on the misalignment angle, 
the hot spot can either be almost fixed or move in a ring-like 



zone around the magnetic pole (e.g. Romanov a et al.|20 03), 
with an angular velocity dependent on the Keplerian veloc- 
ity of the zone of the disk which the matter comes from. For 
the cases in this paper, the misalignment angle is very small 
(0 = 2°) and we observe that the hot spot almost always ro- 
tates (e.g., Fig. [I] right three columns). In Fig. |A1| several test 
cases at various are plotted, showing that the spot move- 
ment is less and less pronounced as increases, as expected. 

This hotspot movement can be visualized in this way: for 
a moment let us imagine the star as a perfectly aligned rotator. 
Matter falling onto the polar caps forms a sort of whirlpool. 
If density enhancements of the flow are formed (see Fig. [I] 
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Figure 3. Position of the hot spots, for the case r c = 1.8 (t* = 5.4ms), a = 0.03, \x = 0.3. This is an example of accretion in the magnetic boundary 
layer regime. The time difference between the files is 0.25P . 




and also Fig. [2] bottom left: the density over the funnel is 
not constant) these density enhancements rotate around the 
magnetic pole with the angular velocity of the zone of the disk 
they come from. 

But when the disk comes closer, either for a weak mag- 
netic field or for a high accretion rate, Rayleigh-Taylor in- 
stabilities appear at the interface between the disk and the 
magnetosphere. These instabilities produce long and narrow 
tongues of matter, that penetrate the magnetosphere and fall 
towards the equator. Matter spilling through the magneto- 
spheric boundaries can, at this point, fall onto the equator or 
be diverted by the magnetic field lines to a random point be- 
tween the equator and the magnetic poles (see Fig. [2] left; see 
also |Kulkarni & Romanova|2008| |2009P . In this case multiple 
hot spots appear, whose duration is normally shorter than in 
the stable case, and whose motion is more complicated, be- 
cause the tongues of matter originating them follow random 
paths before reaching the star. The formation of instabilities 
does not usually prevent the funnel from forming, as Fig. [2] 
shows. The angular velocities of the unstable spots are mea- 
surable from the spot-omega diagram described below and 
thus, as we are going to see in the next sections, they can give 
observable effects on the lightcurves. 

If the disk comes very close to the surface of the star, 
matter can penetrate the field lines forming two ordered an- 
tipodal tongues which fall towards the equator of the star 
(Romanova & Kulkarni 2009), as in Fig. [3] The hot spots pro- 
duced in this way last longer than in the above cases, and 
move with a stable angular velocity corresponding to the Ke- 
plerian velocity at the truncation radius. 



3.2 The spot-omega diagram 

We need a tool to describe the movement of the hot spots 
on the surface of the star in a quantitative way. A tool which 
not only highlights the movements, but also gives us a way to 
measure their angular velocity and their duration. 



Figure 4. The contour levels show the density distribution on the 
star's surface, as seen from the magnetic pole. To draw the spot omega 
diagram, we choose a ring around the magnetic pole, delimited by the 
dashed lines. Then, for every direction cp we plot the average density 
in the area delimited by the ring borders and the directions 4> + dcj) 
and 4> - dcj), obtaining a plot like the bottom of Fig. [5] 



The spot-omega diagram, first introduced in Kulka rni &| 
Romanova ( 2009 ), shows the variation of the azimuthal den- 
sity distribution around the magnetic pole during accretion. 
To obtain it, we choose a ring around the magnetic pole as 
in Fig. [4] whose limits are two angles 1 and 2 with respect 
to the magnetic pole. Then, for every direction 4> around the 
magnetic pole, we calculate the average density of accreting 
matter in the area enclosed by ± A<£ (with A<£ a small step) 
and the borders of the ring. We do this for every time step of 
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our simulation. This gives us the azimuthal density distribu- 
tion as a function of time. 

In the resulting diagram, rotations of the spots appear as 
drifts of the maximum of the flux with time. If we draw a line 
along the drifts, evidently co = Ac/) /At, and so the slopes of 
the curves give us a measure of the angular velocity of the 
moving spots. 

To highlight movements which are longer than one rota- 
tion around the star, we can replicate the plot in the vertical 
direction, going from = to, say, <fi = 8n instead of 2n, 
obtaining for example the plot shown in Fig. [5] 

The spot-omega diagram can be analogously used to 
track the hot spots formed in the equatorial plane, choosing as 
9 1 and 6 2 two angles that enclose the equator, measured with 
respect to the rotational axis instead of the magnetic axis. 



3.3 Measurement of the hot spots' frequency with the 
spot-omega diagram 

In the spot-omega diagram the movement of the spots is 
shown by the variation of position of the flux maxima with 
time. The angular velocity co of the spots is straightforward 
to calculate from the slopes of the curves the moving spots 
produce. The values in the boxes are the angular velocities in 
the frame of reference of a distant observer, in units of co*. 

Before utilizing the frequencies obtained from the spot- 
omega diagram, let us justify its use. Our goal is to show 
that the movements whose frequencies can be observed in 
the spot-omega diagram also produce oscillations in the 
lightcurves, and thus traces in the power spectrum. 

After obtaining the diagram, we independently calcu- 
late the lightcurve produced by the surface of the star. This 
lightcurve is produced in the approximation of perfect black 
body emission of all the kinetic energy of the accreting matter 
and includes corrections for relativistic Doppler effect, light 
bending and time delay ( Kulkarni & Romanova 2005 ). 

We analyze it with a wavelet transform, in order to ver- 
ify if frequencies associated with the spot movement can be 
extracted by means of spectral analysis. 

Here we show an example of stable accretion onto a neu- 
tron star. The parameters chosen for this case are r c = 1.5 
(corresponding to a NS with t* = 4.1ms), = 2°, a — 0.04, 
\x — 2. The results for this case can be summed up as in Fig. [5] 

Not surprisingly, we can see a correlation between the 
slopes observed in the spot-omega diagram and some features 
of the wavelet transform. Moreover, we should take into ac- 
count the fact that the typical time scale of our simulations is a 
few rotations of the star (milliseconds in the case of AMXPs) . 
The wavelet transform only has a very short time to record 
spectral features. In real accreting systems, time scales for av- 
erage density changes in the disk are much longer. Therefore, 
we can reasonably expect to see the emission from the moving 
spots as spectral features, in contrast with the usual assump- 
tions found in the literature. 

As shown in Fig. [6] the velocity of the spot is consistent 
with the Keplerian velocity of the zone of the disk where the 
funnel starts from. 

Thus, so far we have shown that spot motions produce 
visible features in the power spectra, and can be tracked using 
the spot-omega diagram. 



4 MONTE CARLO SIMULATION OF LONG-TERM 
LIGHTCURVES AND PRODUCTION OF QPOS 

In the last section we showed that during accretion, in both 
the stable and the unstable cases, spots can move on the sur- 
face of the star with a different velocity than the star itself. 
What is evident from our simulations is that the moving spots 
are as bright as the ones produced in other simulations, in 
which the spots were fixed. If the fixed spots are able to pro- 
duce the strong and stable signal we observe in Accreting Mil- 
lisecond Pulsars, we believe that the moving ones must be 
visible through particular features in the power spectrum. 

We suggest that QPOs may be among these features. To 
show this, we start by observing that the spot movements 
shown in 3D simulations last a very short time, only a few ro- 
tations of the star, because disk conditions are changing very 
fast. In real accretion disks, disk properties change on much 
longer timescales, say minutes at least. We can safely assume 
that if movements of spots exist, they will probably last much 
longer than in simulations. Nonetheless, we can expect the 
flux of matter not to be constant, and so the spots may be 
sometimes brighter, sometimes less. 

These bright moving spots emit a signal which, from an 
observer point of view, must be modulated with a frequency 
corresponding to the angular velocity of the spots in the ob- 
server's frame. Our way to model this kind of emission is to 
imagine consecutive trains of sinusoids, of variable duration 
and amplitude but constant or quasi-constant frequency. To 
simulate the non-pulsed emission from the disk, we add a 
fixed background. 

Furthermore, we see from simulations that near the v*-\x 
plane, where the infall of matter is slightly more favorable, 
the spots are slightly brighter on average. We model this as an 
additional pulsed signal to add to the rest of the lightcurve. 
This pulsed signal, coming from a fixed point on the star, is 
rotating with the star and will thus have the frequency of the 
star v*. 

At the end of this process, we have a single, continu- 
ous lightcurve in which we can separate three components: 
a fixed background representing the overall average emission 
from the star, a continuous pulsed signal with the frequency 
of the star, and a certain number of wavetrains with the fre- 
quency of the spots. 

As a final icing on the cake, we notice that normally ac- 
creting neutron stars are observed in X-rays, and therefore 
the observed signal is not a continuous lightcurve, but rather 
a series of photons which are detected as "events" by the in- 
struments in X-ray telescopes. To switch from a continuous 
lightcurve to a more realistic series of events, we take advan- 
tage of the fact that common detection rates for this kind of 
objects in real observations are of the order of 10 2 ~ 3 y/s. For 
millisecond pulsars, we have approximately one event per ro- 
tation, and we can thus safely assume these events as rare. 
In fact, these events are very well described by a Poissonian 
variate with an amplitude of probability corresponding to the 
number of photons detected per unit time. 

Then, we obtain a Poissonian random variate with the 
given amplitude and thus create an event series with the same 
statistical properties as the "real" ones. This helps us compare 
the results in our Monte Carlo simulation with the observa- 
tional properties of QPOs from real observations. 
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Figure 5. Frequency analysis for r c = 1.5, a = 0.04, jl = 2, corresponding to a NS with t* = 4.1ms and B = 10 9 G. Top panel: wavelet transform of 
the emission from the surface of the star, as seen by an observer at 45 degrees. Bottom: the spot-omega diagram. The dashed lines lines follow the 
movement around the magnetic pole of the zones at brightest intensity ("hot spots") as explained in j3.2| Values in the boxes show the angular 
velocity in units of angular velocity of the star, and the same values are plotted with the blue dashed lines on the wavelet diagram. The horizontal 
small features in the spot-omega diagram are due to artifacts intrinsic to the cubed-sphere grid for very small misalignment angles. 



4.1 Fixed wavetrains and QPOs 

The movements observed in 3D simulations produce a modu- 
lation in the lightcurve emitted by the star. Let us assume that 
this modulation is purely sinusoidal, but of fixed duration, 
and thus the product of a sinusoid and a rectangular window 
function. It can be shown that the Fourier transform of such a 
sinusoidal wavetrain whose duration is T is a broad peak for 
which Av w 1/T, where Av is the width of the peak. Hence, 
if we say that T is x times the period t of the impulses, then 
the quality factor is: 



V T XT 

Q= — = - = — =x. (10) 

Thus, Eq. [To] effectively says that the Q factor of a broad peak 
in the power spectrum produced by short wavetrains gives us 
an estimate on the duration of the wavetrains themselves. In 
nature, the Q factor of kHz QPOs is such that, if the model we 
propose is correct, the spot movements should last for hun- 
dreds of rotations in some cases. In our 3D simulations the 
hot spots rotate with a stable frequency for a much shorter 
time, but we must consider that the whole simulations only 
lasts tens of rotations, due to numerical problems that do not 
permit to maintain stable disk conditions. In nature the disk 
conditions are, in all probability, much stabler. 



4.2 Description of the model and implementation 

In this section we show in more detail the parameters of our 
model and the procedure to produce the simulated lightcurve. 

4.2. 1 Random number generation 

To generate the random variates we used the very well tested 
functions of the GNU Scientific Library v. 1.10 (Gal assi et al.| 
2004 ) . As random number generator, we chose the luxury ran- 
dom number algorithm of Liischer (ranlxd2), characterized 
by a period of about 10 171 , as implemented in the library itself. 

4.2.2 Parameters 

Every moving spot is represented by a wavetrain described by 
the following parameters: 

• spot frequency v s , corresponding to the angular velocity 
of the spot in the observer frame of reference; its value can be 
constant or vary in a given range; 

• duration AT, whose values can assume a range of values 
between (AT) min and (AT) max ; 

• background B, the amplitude (in y/s), constant, of the 
simulated lightcurve; 

• oscillation amplitude A, constant for a given wavetrain, 
but varying between different wavetrains from A min to A max . 
Its value is expressed as a pulsed fraction of the background. 
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Figure 6. Density and angular velocity distribution around the star at t = 15 (Left) and t = 30 (Right), in the xz plane (Top) and along the x 
direction in the equatorial plane (Bottom). In the plots at the top, blue contour lines represent the omega distribution, while floods represent 
the density distribution. The boxes mark the contour lines corresponding to the velocities seen in Fig. [5] The green continuous line represent the 
magnetic boundary layer (for which matter pressure is equal to the magnetic field pressure) and the dashed line marks the corotation radius. In 
the bottom plot, the black line marks the radius at which matter has the same angular velocities as in the top plots. 



It can also be further modulated as exponentially decreasing 
or Gaussian shaped during a given wavetrain; 

• next spot appearance time t next , corresponding to the 
time between the appearance of two consecutive spots, again 
varying between t next to t next : 

J © min max 3 

• phase 4>, which can assume any value between and 2n. 

Furthermore, we can also choose the frequency of the 
star co*, used to produce the fixed spot emission and also as a 
unit of time to be consistent with the units used in 3D simula- 
tions, and the background level B, constant during the whole 
simulation. 



4.2.3 Procedure 

We first produce the continuous lightcurve by taking the fol- 
lowing steps: 

(i) we choose a frequency v, a sampling time t and a total 
length L of the lightcurve; 

(ii) we choose a pulse shape & (by default sinusoidal) and 
a window function IV (typically flat, Gaussian or exponen- 
tially decaying as in Fig. [7} giving the overall intensity vari- 
ation of the wave trains. It is of particular importance to also 
choose if we want a constant average intensity of the source 
(and so J &{t)dt = 0) or if we want the spots to create an 
additional luminosity ("burst-like"); 



Train shapes 

1 j iii 

— Constant 

— Exponential 

— Gaussian 



700 - 




0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 



Figure 7. Examples of possible synthetic wavetrains: we can notice 
the constant sinusoidal modulation with the frequency of the star v* 
and several different behaviors of the intensity of the hot spot, with 
its modulation at the frequency of the spots v s . 



(iii) we generate a random value (distributed according to 
a uniform or Gaussian variate) for A and AT; 

(iv) we generate a random value (distributed according to 
a uniform, Gaussian or Exponential decay-like variate) for 

j-next . 

L 5 
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(v) we produce a random value for 4> ', 

(vi) we generate a wave train of duration AT from t n 
according to the formula: 



a(t)=A&>(t + <j)T s )W(t) 



(11) 



(vii) we go back to |(iii)| if t < L 

(viii) to the lightcurve obtained, we add a constant back- 
ground and a continuous sinusoidal modulation with fre- 
quency CO*. 

At this point, we have the complete simulated lightcurve as 
it would be seen by an optical telescope. We only need the 
following few steps to obtain a series of events as in real X- 
ray observations: 

(i) for every t we produce a random integer number N events , 
distributed according to a Poissonian variate with amplitude 
TA{t), corresponding to the expected number of events in the 
time t; 

(ii) we write in the output file the number of events N events 
"recorded" at the time t. 



4.3 Spectral analysis 

Fourier Transform is the main technique used to find period- 
icities in time series. The Fourier transform, in fact, gives a 
representation of the input function in term of its harmonic 
content. With Fourier transforms it is possible to find period- 
icities in a time series even if the amplitude of the periodic 
component is very small. 

We performed a Fast Fo urier Transferred on the data by 
means of the FFTW3 library ( [Frigo & Johnson|2005] ), and ob- 
tained the power spectrum according to the formula 



P J = 



x 2 + y 2 
N 



(12) 



where x and y are the real and imaginary parts, respectively, 
of the Fourier transform of the time series, and N is the total 
number of simulated events. This is the standard normaliza- 
tion proposed by Leah y et al] ( |1983| ), which is particularly 
useful in the case of Poissonian events. In fact, using this nor- 
malization the power spectrum of a constant background fol- 
lows a x 2 distribution with 2N bin degrees of freedom, where 
N bin is the total number of bins in the spectrum. In this way 
we have a very solid statistical framework to help us interpret 
our results. 

In order to improve the signal to noise ratio, we per- 
formed the FFT on M consecutive chunks of data and obtained 
a single FFT from the average of the chunks. With the cho- 
sen normalization, the backg round has an avera ge of 2 and a 
standard deviation of 2/Vm ( [van der Klis|l988| ). 



4.4 Example of Monte Carlo-produced QPO 

The case in Fig. [8] is corresponding to a spin period = 3ms 
and a hot spot period t 5 = 2ms; we used a duration of the spot 
motion of AT = 50ms and t next — 100ms. The pulsed fraction 



1 There are several algorithms to calculate Fourier transforms of sam- 
pled data. The fastest ones, which exploit the properties of Fourier 
transforms and also of the data representation in computers, are 
known as Fast Fourier Transforms. 



-e-47i 




Figure 11. Spot-omega diagram for the case r c = 1.8 (T* = 5.4ms), 
a = 0.03, Jx = 0.9. In this case is possible to notice the coexistence of 
two motions of the hot spots, with two different scales of velocity. 



of the moving spots was 20%, while that of the fixed spot was 
10%. We can see the QPO peak at 1.5 and the sharp peak at 
the frequency of the star. 

On the right hand side of the same figure, the QPO peak 
is fitted and some details are provided: as one can see, the 
peak is well fitted by a Gaussian with cr = O.Olv*. 



5 CORRELATIONS WITH THE MASS ACCRETION RATE 

In §[3] and [4] we showed that movements of the hot spots can 
produce QPOs. In this section we show that those movements 
are not a fortuitous event, but in fact they present a pre- 
dictable behavior. To show this, we study the change on the 
angular velocity of the hot spots as the mass accretion rate 
changes, performing a number of new simulations with differ- 
ent values of jl using Eq.[8]to find the corresponding values of 
M. We then represent the frequencies observed from the spot- 
omega diagrams against the corresponding values of jl and M 
in plots like those in Figs[9]and[To] To be complete, we include 
in the study all the frequencies visible in the spot-omega dia- 
grams for all cases. We do not perform any choices on the fre- 
quencies to plot. Multiple frequencies appear in some cases 
simultaneously, in others at different times. The size of the 
bullets in the plots is proportional to the time the plotted fea- 
tures appear at during the runs. In this way is possible to see 
whether multiple features appear together or not. 
In Figs.[9]and : 



10 



) the Vqp s obtained from the spot-omega 
diagrams are plotted against the mass accretion rate and 
the dimensionless parameter fl (which determines the size 
of the magnetosphere) for several cases. As previously said, 
we choose to model a neutron star with B = 10 9 G, R = 10km, 
M = 1.4M . According to Eq.[8]we can use p, to change M , in- 
vestigating the variation of v QP0 with respect to the measured 
value of M. 

The case in Fig. [9] describes an NS with t* = 4.1ms. This 
is a case of almost stable accretion, except for high values 
of the accretion rate. As the accretion rate grows, the mag- 
netosphere becomes smaller and smaller, and the hot spots 
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Power spectrum 

T 




Power spectrum, detail 




Figure 8. (.Left) Spectrum of the simulated time series, with t* = 3ms, t 5 = 2ms, AT = 50ms, t n 



100ms, A = 20%, B = 300y /s. In this plot 



we can observe the QPO at v = 1.5V* and the peak at the frequency of the star (v = 1). The length of the simulated time series was 10000s and 
the spectrum was obtained averaging the spectra of 16384 subintervals. (Right) Gaussian fit of the simulated QPO 



700 



vqpo vs.log 10 M (z/* = 245.70Hz) 



-10.0 
600 

500 - 

M 400 
O 

^ 300 - 
200 
100 - 





-9.0 



~l 

-8.5 



-8.0 



A 
▲ 



A 
A 



log^MiMey- 1 ) 



vqpo vs. Jjl (is* = 245.70Hz) 




0.5 1.0 1.5 



Figure 9. Left: correlation between the frequencies vq P0 observed in the spot-omega diagram and the mass accretion rate for the case r c = 1.5 
(T*=4.1ms), a = 0.03. The scale at the top is for a neutron star with B = 5 • 10 8 G, the one at the bottom for B = 10 8 G. Right: correlation with 
the dimensionless parameter Jx in the same case. The dashed line marks the frequency of the star. The colors indicate fx, in order to compare the 
position of the points in the two plots, and the size of the mark is proportional to the time of appearance of the features in the simulation. 



move around the magnetic pole with an increasing frequency, 
as expected. When M reaches higher values the accretion is 
not stable anymore, and multiple hot spots appear, originat- 
ing the multiple frequencies observed in Fig. [9] In this case all 
the simulations were short compared to the ones in the next 
paragraph because of the appearance of numerical instabili- 
ties, and so all frequencies are observed at an early time. 

The case in Fig. [To] describes unstable behavior, like we 
can see in Fig. [2j In this case we can observe a peculiar be- 
havior: the hot spots in the ring around the magnetic pole and 
the ones produced by instabilities move with different angu- 
lar velocities. This movement is visible also from the spot- 
omega diagram in Fig. [Tl] This produces two traces in the 
spot-omega diagram, whose frequencies present the same ap- 
proximate behavior, as M changes, as the ones seen in the 
first case. The frequency separation of the two peaks is in the 
correct range of values (200 - 300Hz) to suggest that the twin 
kHz QPO peaks often observed in nature could be originated 
in a similar way. A confirmation of this idea can be the fact 



that the upper kHz QPO is usually much less coherent than 
the lower one fldi Salvo et~aLl|2003[ |Barret et al1|2006| in- 
vestigate this particular aspect for a wide sample of LMXBs) . 
In our simulations, the "upper" peak is produced by instabili- 
ties while the "lower" one is produced by the motion around 
the magnetic pole. Normally, unstable accretion is character- 
ized by a much smaller coherence than the stable drift around 
the magnetic pole, except in the strongly unstable regime in 
which the two antipodal tongues maintain a constant angu- 
lar velocity for a relatively long time (Romanova & Kulkarni 
[2009] ). 

Both frequencies relate to the zone near the inner disk. 
The instabilities form at the edge of the disk, and hotspots 
produced by them have higher frequencies than those of the 
hotspots in the polar region, which are produced by mat- 
ter captured by the field lines at a slightly larger distance 
(see Fig. [6] Fig. [2}. The two frequencies move together as 
M changes because both of them depend on the distance of 
the disk from the star. 
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Figure 10. v QP0 variation with respect to log 10 M and /I for r c = 1.8 (T*=5.4ms). As before, the scale at the top of the u Vq P0 vs. M" plot is for a 
neutron star with B = 5 • 10 8 G, the one at the bottom for B = 10 8 G. Top and center cases have a = 0.02 and a = 0.03 respectively, while at the 
bottom the two cases are plotted simultaneously. 



6 DISCUSSION 

We developed a model of quasi-periodic oscillations based on 
the spot movements on the surface of the accreting neutron 
star with a slightly misaligned dipole magnetic field. 3D MHD 
simulations have shown that spots may rotate faster/slower 
than the star in cases of stable or unstable accretion (see also 
Rom anova et al,||2003| |2004[ |Kulkarni & Romanova||2008l 
2009) and in the magnetic boundary layer regime ( Romanova 



& Kulkarni 2009 ). We have shown that in all above cases, 
the spot movements develop ordered oscillations in the light- 
curves. Modeling these oscillations as short-lived sinusoidal 
trains, we showed that if they are powerful enough to be seen, 
their fingerprint is a QPO. Given the range of values of the fre- 
quencies involved, we hypothesize that kHz QPOs could be 
produced in this way. 

There are several features of QPOs that can give con- 
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firmation to this model. First of all, the general trend with 
the mass accretion rate, as shown in ^5] is as it is expected. 
Though it is not generally true that luminosity and mass ac- 
cretion rate are proportional to each other, we can assume it 
to be true on a given source and for a short timescale. This is 
the case of the simulations in this paper, where only one pa- 
rameter (the mass accretion rate) changed between the vari- 
ous simulations. 

The rms amplitude of real life QPOs has a strong depen- 
dence on energy (Ber ger et al.||1996] ), suggesting that the 
mechanism of emission is different from that of the back- 
ground. In our opinion, this is because the QPOs are produced 
on the surface, or perhaps just above the surface by means of 
shocks in the accretion column similarly to what happens in 
AMXPs (see for example Gierrins kTet al.|20 02 ), while most of 
the background is produced by the disk. The rms amplitude 
is also generally higher in atoll sources (up to 20%) than in Z 
sources (2-5%) ( [Jonker et al.|2001| ). The reason for this dif- 
ference is probably related to the higher accretion rate of Z 
sources. 

In some cases 3D simulations show that two frequencies 
appear at the same time. This numerical discovery is impor- 
tant for possible interpretation of two QPO peaks observed 
in a number of LMXBs. In fact, this mechanism for the pro- 
duction of multiple QPOs is able to explain more than the 
mere appearance of two peaks in the Power Spectrum. When 
double QPOs appear in observations, they normally present 
different behavior as regards their Q factor (di S alvo et al] 
2003] |Barret et"aL||2006| ), and hence (according to our as- 
sumption) their coherence time with the upper peak gener- 
ally less coherent; in our model, this is due to the fact that the 
lower one comes from a smooth motion around the magnetic 
pole, while the upper one is produced by instability tongues, 
shorter-lasting and thus with a lower coherence. The term 
"unstable" should not mislead the reader though: the unsta- 
ble regime is very common in simulations, expecially for fast 
rotators, also at small values of the mass accretion rate, as 
shown in Fig. [To] A higher coherence of the upper kHz QPO, 
in our model, could in principle be produced by a very high 
accretion rate and the onset of the magnetic boundary layer 
regime, in which instabilities produce long-lasting opposite 
tongues in the equatorial plane. During this regime though, it 
is very probable that the zone around the surface of the neu- 
tron star is optically thick. Just as observations show, multiple 
QPOs can appear simultaneously or at different times. 

As regards the models predicting the upper limit of v u as 
the Keplerian frequency at the ISCO ( [Miller et aT]|1998| ), it 
would be interesting to investigate the matter further. For the 
values of mass, magnetic field and radius used in the present 
job, the ISCO is never reached by the inner disk during QPO 
production with the two channels described because of the 
onset of the magnetic boundary layer regime. Further inves- 
tigation should include the use of more extreme values of the 
mass in order to study the behavior of the inner disk when 
the ISCO is considerably far from the surface. 

Another interesting point is the v - L x correlation ob- 
served in many papers (for example Yu et al. 1997; Mendez 
|et al.|1998] ). The correlation only holds for a given source and 
on time scales shorter than a few hours. The main reason why 
it does not hold across sources can easily be found in the mag- 
netic field: for a given value of the mass accretion rate, differ- 
ent values of the magnetic field correspond to a different size 



of the magnetosphere and, as a consequence, to different val- 
ues of the frequency of the moving spots. Again, this is evident 
from our simulations: the exact same frequency in figure [9] is 
obtained for M w 4 • \0~ ll M Q y~ l in a system with B = 10 8 G 
and for M = lO~ 9 M y _1 in a system with B = 5 • 10 8 G. But the 
simulations presented here show only one of the many pos- 
sible configurations of the disk around the star. For a given 
source, it is possible that changes in the structure (density 
distribution for example) of the disk produce a different be- 
havior of the moving hot spots, even with the same mass ac- 
cretion rate. Besides, the relationship between luminosity and 
mass accretion rate is not clear. This point will be investigated 
more deeply in future papers, by simulating systems similar to 
the one shown here but with different initial disk conditions. 

The results presented here are based on simulations of 
neutron stars with very small misalignment angles. An equiv- 
alent discussion can be made at slightly larger misalignment 
angles, where the behavior of the spots does not change con- 
siderably (see Fig. |A1[ but also |Kulkarni & Ro manova 2008 , 
2009; Romanova et al. 2008]), and obtain almost the same 
features in the power spectrum. In this case though, from sim- 
ple geometrical considerations it is reasonable to expect some 
kind of modulation at the frequency of the star. The least we 
can expect is a broad Lorentzian-like feature at the frequency 
of the star ([Burderi et al.|1993[[Lazzati & Stella|1997[|Menna| 



et al. 2002 ) just from the action of the rotating spot in the 
polar region, due to the rotation of the star. Its quality fac- 
tor Q would be comparable to that of the moving hotspots 
(Q*/Q sp ot w vJVf according to the considerations of 



The observation of such a feature, even very dim, in observa- 
tions where the lower kHz QPO (the one from the funnel flow 
in this model) is present, would give interesting evidence in 
favor of this model. 

As a final remark, we notice that the difference between 
the frequencies of the two QPOs is 200 - 300Hz, close to the 
frequency of the star. It is generally known that for relatively 
slow rotators the frequency difference of QPOs is similar to 
the frequency of the star. Until about two years ago, it seemed 
that the Av was approximately equal to v* for slow rotators 
and 0.5v* for fast rotators. jMendez & Belloni| ( |2007| ), ana- 
lyzing the matter in detail, find that this distinction is not 
so clear, and that data can be interpreted and fitted equally 
well by values of Av around 300Hz for all sources, with the 
caveat ( [van Straaten et al."l20 05') of a few AMXPs showing all 
variability components (and thus also Av) shifted down by a 
factor of ~ 1.5. The range of values of Av found in this work 
are compatible with both hypotheses, because the simulated 
system is a slow rotator and Av is both near the frequency 
of the star and 300Hz. We are planning to investigate the 
matter with simulations of faster rotators. Nevertheless, this 
mechanism for the production of double QPOs does not need 
the two frequencies to be origined from one another, giving a 
fixed value of 5v, which was the main caveat of a somewhat 
similar model, the sonic point beat frequency model by Miller 
letaLl ( [19981 ). 

From what we observe in this work, the two QPOs are 
two separate effects of the interaction between the rotating 
magnetosphere and the disk. The funnel flow comes from a 
zone just behind the inner radius where the magnetic field 
is strong enough to capture matter and make it drift towards 
the pole. The instabilities form right at the inner radius, where 
magnetic field energy density and ram pressure are equal (the 
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green line in Fig. [TJ . The results of this paper can give insight 
into the matter: the frequencies found in observations are a 
probe of the interaction between the magnetic field and the 
disk. 
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APPENDIX A: ON THE BEHAVIOR OF THE SPOTS AT 
DIFFERENT MISALIGNMENT ANGLES 

In Fig. |Al] we can compare the behavior of the hot spots at 
different misalignment angles. As can be seen in the picture, 
the spot movement is not limited to the case with very small 
misalignment angles, being present even for = 15°. 
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Figure Al. Hot spot shape for r c = 1.8 (t* = 5.4ms), a = 0.03, jl = 0.5, and (Top) 9 = 5°, (Middle) 9 = 15°, (Bottom) 9 =30°. For small 
misalignment angles, the hot spot moves with an approximately constant intensity around the magnetic pole. For bigger angles, the hot spot 
tends to stay fixed or to vary its intensity from a maximum value in the preferred position (which depends on the misalignment angle and on M) 
to a minimum value at the opposite side with respect to the magnetic pole. 
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